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Abstract 

Numerical simulations of self-gravitating systems are generally based on N-body codes, which 
solve the equations of motion of a large number of interacting particles. This approach suffers 
from poor statistical sampling in regions of low density. In contrast, Vlasov codes, by meshing 
the entire phase space, can reach higher accuracy irrespective of the density. Here, we performed 
one-dimensional Vlasov simulations of a long-standing cosmological problem, namely the fractal 
properties of an expanding Einstein-de Sitter universe in Newtonian gravity. The N-body results 
were conhrmed for high-density regions and extended to regions of low matter density, where the 
N-body approach usually fails. 
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I. INTRODUCTION 


In the present epoch, the observable universe is extremely inhomogeneous. Taking galaxies 
for the essential elements, we see them grouped in clusters that are, in turn, grouped in 
superclusters, and interlaced with enormous voids. The recent detection of the local su¬ 
percluster Laniakea [Tj exemplihes this scenario. Consideration of these gross features lead 
Mandelbrot [2] and Pietronero [3], among others, to conjecture that the universe is a fractal, 
at least at some intermediate scales. An early thinker along these lines was de Vaucouleurs 
[1] . Since cosmological theory demands that the universe is homogenous at sufficiently large 
scales, the search for the transition to homogeneity has been a focus of recent investigations 
W- Partial evidence for the fractal geometry is provided by the two-body correlation 
function computed from recent large-scale galaxy surveys like the Sloan or 2-degree m, 
which exhibits power-law behavior over a hnite range of scales. However, it is difficult to 
determine the mechanism and evolution of such scaling behavior from either observation or 
three-dimensional (3D) N-body simulations. To obtain a more complete picture, investiga¬ 
tors have turned to lower-dimensional models where a more precise representation of the 
gravitational held is possible over all scales. 

A great deal of work on structure formation in the universe has been accomplished using 
Newtonian ID models (for a review see [TT], and for more recent work [T^HTTj i. The link 
between ID and 3D cosmology models was discussed in [TS]. Nevertheless, N-body simu¬ 
lations, even in ID, suffer from an intrinsic undersampling of the phase space because of 
the hnite number of particles used in the codes [16]. As in reality the number of bodies 
is virtually inhnite, one should instead solve the mean-held Vlasov equation (which is the 
V —)■ oo limit of the N-body model) for the phase-space distribution coupled to the Poisson 
equation for the gravitational held. This is a very demanding computational task, both in 
terms of run duration and memory storage, particularly for situations where high accuracy 
is necessary to resolve the intricate phase space structures that develop over time. However, 
present computers now make this approach feasible, if not in 3D at least for ID models. 

This work is devoted to the presentation of the hrst cosmological results obtained with a 
ID Vlasov approach. The paper is organized as follows. In Sec. [IT] we will introduce a 
set of scaled variables (comoving coordinates) that are particularly adapted to the Vlasov 
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approach. Some details of the numerical algorithm used to solve the Vlasov-Poisson equa¬ 


tions are provided in Sec. Ill The numerical results are presented and analyzed in Sec. IV 
General conclusions are presented in Sec. |Vj 


II. MODEL AND SCALED VARIABLES 

Let us consider a highly symmetric expanding distribution of matter. Its gravitational held 
has only one component Er{r,t) which depends on time and on a single spatial variable 
r. The symmetry could be, for instance, spherical or planar. In the planar case, originally 
developed by Rouet and Feix (RF) [T^ EH] and later expanded by Miller and Rouet [T9H2T], 
the system is constituted of many parallel expanding planar sheets whose surface density 
decreases following the expansion law. For spherical symmetry (the so-called quintic model 
|22j ) the system is composed of concentric spherical shells. 

The equation of motion of a particle in such a Newtonian gravitational held reads as 

§ = M ( 1 ) 

where r{t) is a spatial position in the expanding universe. In the mean-held limit, the 
gravitational held is a solution of the Poisson equation 

Vr ■ E = —AirGp, (2) 

where p{r, t) is the matter density. In order to account for the expansion, we transform 
space and time according to: 


r = cm, (3) 

dt = A^{t)d6, (4) 

where ^ is a comoving spatial coordinate. With this scaling, the velocity variable is trans¬ 
formed as 

V = CAm + Ci, (5) 

where v = dr/dt and r] = d^/dd (a dot denotes diherentiation with respect to the time 
t). Here A{t) and C{t) are two strictly positive functions of time. The general equation of 
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motion in the scaled variables reads as 


de^ 


+ 2A^ 




dd 


.4^ 


As 

C3 ’ 


where S{^,9) = C‘^(t)Er{r,t) is the scaled gravitational held, satisfying 


Vg ■ £ = -dvrGp, 

and p{^,0) = C^(t)p{r,t) so that the total mass is preserved. 
For the scaling functions, we use the following forms: 

.4^(4) = ; C(i) = (t/t„r' 


( 6 ) 

( 7 ) 

( 8 ) 


A. Standard scaling 


The standard scaling [13 [H] uses (3 = 1 and 7 = 2/3. With this choice, all coefficients in 
Eq. ([^ become time-independent: 


do^ m„d0 9tf 


( 9 ) 


For a constant density p = po, Poisson’s equation ([^ can be solved exactly in a d-dimensional 
space to give the gravitational held £ = —AnGpo^/d = —coj^/d, where uj = y/AnGpo is 
the Jeans frequency. At equilibrium, the gravitational held must exactly cancel the inverse 
harmonic term on the left-hand side of Eq. 0 . This provides the relationship between to 
and Uj: 


Ujto = -d. 


Therefore, Eq. (1^ can be written as: 


A^{t) = {aujjtY ; G(t) = {aujty, 


( 10 ) 


( 11 ) 


with a = 1 /( 0 ;jto) = 3/-\/M, so that the scaled and unsealed co-ordinates coincide at t = to- 
Then, Eq. ([^ becomes 

( 12 ) 


, UJJ di _ ^ ^ 


d6‘^ d9 d 

The RF model (considered here) has planar symmetry and is therefore essentially one 
dimensional {d = 1). The corresponding scaled Poisson equation is also ID: d^S = —AnGp 
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For the quintic model [22], which corresponds to a spherically-symmetric expanding universe, 
we have d = 3. If we consider a planar perturbation in the scaled universe, we can still use 
the ID Poisson equation as above; however, the factor in front of the background term (third 
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term on the left-hand side) of Eq. (12) must be modihed as follows: —)■ cjj, in order to 
allow for a steady state at equilibrium. (Note that there is no such change in the RF model, 
because both the original system and the perturbation have planar symmetry). With this 
substitution, the scaled equations of motion of the RF and quintic models only differ in the 
coefficient of the friction term, and can be written as 

, 1 de 


dd2 y^dd 




(13) 


In the above equation we also introduced nondimensional variables, whereby the scaled time 
6 is normalized to the inverse Jeans frequency, the scaled space coordinate ^ is normalized to 
an arbitrary length A, and the scaled gravitational held S to Xuj. We keep the same symbols 
for the nondimensional variables, which will be used throughout the rest of this work. The 
nondimensional Poisson equation reads as: d^S = —p where the density is normalized to po- 


According to Eq. (13), if the universe is homogeneous and strictly follows the expansion 
factor C{t), then it will be static in the scaled variables, with a constant (nondimensional) 
density equal to unity and corresponding gravitational held S = —However, this is an un¬ 
stable equilibrium which, when slightly perturbed, evolves towards a highly inhomogeneous 
distribution of matter with complex features. 


B. New scaling 

It is important to note that, of the two exponents fd and 7 in Eq. (|^, only 7 has some 
physical bearing: it represents the rate of expansion of a self-similar Einstein-de Sitter 
universe. Instead, the exponent (d = 1 was chosen on purely utilitarian grounds to render 
the scaled equations autonomous. This is an appropriate choice for N-body simulations, 
because the ID equations of motion can be solved exactly to machine precision, but a poor 
choice for a grid-based Vlasov code. Indeed, using this scaling, the transformed velocity rj 
grows exponentially in time, as was shown by N-body simulations (see Appendix]^. This 
is a serious problem for Vlasov simulations, which solve the Vlasov equation on a hxed grid 
that covers the entire relevant phase space. If the velocities keep growing, one would need 
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to mesh an increasingly large phase space in the (^, t]) variables, soon reaching memory and 
compntation time limits. 


The important point is that we can choose a different valne of the exponent f3 (and thns 
rescale time and velocity in a different way) so that the scaled phase space stays bonnded for 
the entire dnration of the rnn. This can be achieved by choosing: (3 = (So; — 1)/(3 q;) ~ 0.84 
in Eq. (|^ (details of the calcnlations are given in the Appendix]^. With this valne, and 


nsing the same normalization as in Eq. (|13|), one obtains the scaled eqnation of motion: 


K di e 


(14) 


do^ /i(0) de /i2(0) ’ 

where K = {3 + a/ 2)/(3\/2), and fi{6) = 6/3 + 1. Note that the scaled time 9 depends on 
the valne of /?, which is not the same for Eq. ([I^ (standard scaling, (3 = 1) and for Eq. 


(14) (new scaling, (3 ~ 0.84). Therefore, the two scaled times are not the same and their 
relationship is given in the Appendix 


The Vlasov-Poisson eqnations corresponding to Eq. (|14|) reads as follows: 

dF 


dF 


d 


E 


T] 




F- 


K 


r]F ] =0, 


dr] y]F{9) ]i{9) 

0 

F{^,V,9)dr], 

O 

where F{^,ri,9) is the distribntion fnnction in the rescaled phase space. Note that, by 


^ = 1 


( 15 ) 

( 16 ) 


defining £ = £ + ^, the harmonic term in Eq. (14) has been incorporated into Poisson’s 


eqnation (16). 


III. NUMERICAL METHOD AND INITIAL CONDITIONS 


We solved nnmerically the set of eqnations (15)-(16) with periodic boundary conditions in 


the scaled spatial variable Vlasov codes work by covering the entire phase space {^, 1 ]) 
with a uniformly spaced grid. The distribution function F is pushed in time using a split- 
operator scheme that treats the space and the velocity coordinates separately [23] . The time 
integration between 9 and 9 -I- A9 is performed in three steps. First we solve the equation 


dF 



0 


( 17 ) 
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whose exact solution is just a rigid shift of rjAO in position space: 9 + A9) = 


riA9,ri,9). Then, the gravitational held S is obtained through Poisson’s equation (16). 
Finally, we solve the equation 


dF 


d 


£ 


K 


r]F = 0. 


( 18 ) 


dr] y]j?{6) ]i{6) 

Because of the presence of the friction term and the time-dependent coefficients, this step 


is not standard. However, Eq. (18) can also be solved exactly (considering £ constant), by 
integrating the characteristic: 

dr] £ 

M ^ ]F{e) ~ ~]i{e) ‘ 

which has the following solution 




(19) 


Ar] = r]{e + AO) - r]{e) = C,vi0) + Ce£, 


where Cr, = {B - 1), 


and 


Ce 



— 1 ----— 

FT-1/3 1 + 0/3’ 


B = 


l + (0 + A0)/3 
1 + 0/3 


The solution of Eq. (18) is then: F{^,r],9 + A9) = F{^,r] — Ar],9). 


Interpolations on the phase-space grid are performed using an accurate hnite-volume algo¬ 
rithm m that preserves the positivity of the distribution function. Vlasov codes display 
very low numerical noise even in regions where the matter density is rarehed, which is 
where N-body codes would be most noisy because of poor statistical sampling. The Vlasov 
approach is widely employed in plasma physics, and has been used occasionally to study 
self-gravitating systems |25l [26] , but was never applied to cosmological simulations. A com¬ 
plementary approach to either N-body or Vlasov codes is provided by the water-bag method 


Equations (15)-(16) were solved for 0 < .^ < L and —1 
Vmax = 15. For the present results, we used = 2^^ 
= 1000 points in velocity space. 


max <9 < Vmax, with L = lO^TT and 
points in the spatial coordinate and 


The initial condition is a “cold” Maxwellian in velocity space, with variance {r]'^) = 0.01. 
The initial matter density p{^) = f Fdi] is so constructed as to display a power-law spectrum 


7 














of the type: P{k) = ~ k^, where k is the wave number in ^-space. Initial power spectra 

of the form P{k) ~ /c”, with n G [0,4], were used in a number of earlier works on structure 
formation 

The case P{k) ~ k^ produces the ID version of the Harrison-Zeldovich spectrum from 
the assumption of scale-free potential fluctuations [2H1 122] • Following inflation, density 
fluctuations in the universe can be modeled as a Gaussian random field with a scale-free 
(power law) power spectrum P{k). In a 3D universe, the exponent corresponding to a scale- 
free potential is unity, i.e., P{k) oc k. To see this, we expand the gravitational potential <I> 
in a Fourier series [28l |29] : 


^{f,t) = V ^ <l>^exp(ifc ■ r^, (20) 

k 

where V is the volume. Then it can be shown via the Poisson equation that 

^ ^ f V. ^ f ^ f ^.(.n .). (21) 

Consequently, if P(k) oc k, the potential fluctuations are scale-invariant on a logarithmic 
scale. Initial conditions for 3D simulations of the expanding universe are guided by these 
considerations. 


Similarly, in ID we have: 

^ ^ f ( 22 ) 

Requiring scale-free fluctuations for the potential then yields P{k) oc k^, which is the initial 
spectrum that we took for our ID simulations. 


In order to accelerate the early evolution, this initial condition was hrst allowed to evolve 


for a short time according to Eqs. (15)-(16) where the scale factor /i(6*) was set equal to 


unity. Once the fluctuations have reached a sufficiently high level, the correct scaling was 
applied and the initial time was reset to 6* = 0. 
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FIG. 1: Phase-space distribution functions (left) and matter densities (right) at time 6 = 325. 
The top panels show the entire domain, the center and bottom panels show consecutive zooms. 
The contour levels of the distribution functions are distributed logarithmically in the interval: 
10“® < F < 1. The actual maximum value of F is around Fmax = 185. 

IV. RESULTS 

A. Phase space and matter density 

Figure shows the phase-space distribution function and corresponding matter densities p 
(normalized to unity) at a later time. As expected the velocity domain remains bounded, 
so that the x points that mesh the phase space are used in an optimized way. The 
distribution function clearly displays a hierarchical structure at different scales, with small 
clusters orbiting each other to form larger clusters, which in turn also revolve around each 
other. This hierarchy is at the basis of the fractal structure observed with N-body simulations 
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and discussed later in this work. The density displays many spikes, which become narrower 
and higher as time elapses. These spikes are even more apparent on a semi-log plot of the 
density (Fig. [^. This structure is similar to that observed with N-body codes. 


Note that the straight segments in the phase-space plots of Fig. (see top left panel) all 
have approximately the same positive slope, and correspond to regions of low matter density 
(“voids”). Similar behavior was seen in N-body simulations and it represents regions that 
are devoid of particle crossings [I9H2I]. The slope of these segments can be estimated using 


Eq. (14), where we neglect the gravitational held £ because in the relevant regions the 
density is low. We seek for a solution of the type ^{6) = (1 -|- 9/3)'^. Substituting this 


expression into Eq. (14), we hnd that the parameter 7 must satisfy the algebraic equation: 


7^ + {3K - 1)7 - 9 = 0, 


(23) 


where 3K — 1 = 3/a/ 2. Equation (23) has the positive root 7 = 3/a/ 2 (the other root is 
negative and the corresponding solution is quickly damped away). Then the ratio between 
the phase space variables rj = d^/d6 and ^ becomes: 


7 ^ 7 

e 3 + 0- 


(24) 


For 6 = 325, we obtain a slope t]/^ ^ 0.0065. This is very close to the slope observed in Fig. 

as can be deduced for instance from the segment on the left of the top left panel. We also 
checked that, for large times, the observed slope decreases as 1/6, in accordance with Eq. 

m. 


B. Power spectrum 

The development of hierarchical (scaling) behavior can be tracked by the evolution of the 
power spectra. Here we show in Fig. [^the power spectrum of the matter density as 
a function of the wave number kj = {271/L)j G [2 x 10“^, 3.28], with j = 1... N/2. Rather 
quickly, a decreasing power-law spectrum builds up, with a slope roughly equal to —0.53, 
not far from the value —0.45 observed for N-body simulations of the RF model with the 
same initial spectrum [IT]. The range of the power-law region {kmin^kmax) increases with 
time, with kmin getting smaller and smaller while k^ax remains roughly constant. The steep 
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FIG. 2: Matter densities at time 0 = 325 on a semi-logarithmic scale. 

decrease at A; > kmax is due to numerical diffusion. The observed slope is also consistent with 
recent predictions [121 [S] which, when applied to our simulations, yield a slope of —0.57 
|37j . Benhaiem et al. [T5] found that this result also holds for 3D cosmology in scale-free 
models. 

The power spectrum is a useful indicator of the difference between the Vlasov and N-body 
results. As already noted, the Vlasov power spectrum (Fig. has a slope similar to 
that observed in the N-body case. To gain further insight, we re-analyze the spectrum by 
performing different cut-offs, either removing the low or high values of the matter density. 

Let us hrst consider the high-density power spectrum. In Fig. we show the spectrum 
obtained by considering only the values of the density that are above a certain pmin (values 
that are below this threshold are removed). We observe that the slope becomes less steep 
with increasing cut-off, and is almost flat for p^in = 20. Of course, this is a huge cut-off, 
and we chose to show these values only to point out the general trend. Nevertheless, this 
observation may explain why the observed N-body spectrum is slightly less steep (slope ~ 
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FIG. 3: (Color online). Power spectrum of the matter density \pk\‘^ for different times from 6 = 0 
to 9 = 300. Later times correspond to larger values of \pk\‘^. A moving average over 41 points is 
taken in order to smooth the fluctuations. 

—0.45) than the corresponding fnll Vlasov resnlt (slope ~ —0.53): the N-body spectrnm 
lacks the contribntion from the low-density regions, which tend to steepen the spectrnm as 
seen in Fig. Consistently with this reasoning, the Vlasov spectrnm (which inclndes both 
high and low densities) is closer to the analytical estimate of Ref. [13] (slope ~ —0.57). 
It is also interesting to estimate the valne of the cnt-off that yields a slope similar to that 
observed in the N-body resnlts, i.e., —0.45. We fonnd that the reqnired cnt-off is ~ 4. 
Notice that this is a relatively small threshold in onr nnits, as 94.5% of the matter density 
is above that valne. 

Conversely, if we consider only the low valnes of the density, the observed spectrnm is 
signihcantly steeper, as is shown in Fig. All in all, it appears that the total spectrnm 
resnlts from the combination of a steeper (for low densities) and a flatter (for high densities) 
cnrve. N-body codes only captnre high-density regions, and therefore yield a spectrnm with 
a slope slightly smaller than the theoretical estimate. 

These resnlts may be nsefnl, for instance, to improve onr nnderstanding of the geometry 
and distribntion of particles in voids, which is a cosmological problem of cnrrent interest 
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FIG. 4: High-density power spectrum at 0 = 300. For each curve, density values below the 
corresponding Pmin have been removed before computing the spectrum. 


C. Fractal dimension 

The clustering of the phase-space (Fig. and the power law observed in the density 
spectrum (Fig. point to an underlying fractal structure of the matter distribution, as was 
the case for the N-body simulations |T9]. Box counting is the method most frequently used 
to analyze the properties of a fractal structure [32]. Here, this method is used to determine 
the generalized fractal dimension Dg in ,^-space, also known as the Renyi dimension. The 
system (0 < ^ < L) is covered with boxes of length i of decreasing size: i = i = T/4, 
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FIG. 5: Low-density power spectrum at 0 = 300. For each curve, density values above the corre¬ 
sponding praax have been removed before computing the spectrum. 


and so on. The fractal dimension is defined as 


D, 


lim 


log(^.m,^; 
g — 1 ^^*0 log(£) 
Y,i mi log(mi 


, for g 7 ^ 1 


lim 

£-s>0 


log(£) 


for g = 1, 


( 25 ) 

( 26 ) 


where mi{i) = p{C)d^lmtot represents the proportion of mass contained in the i-th 

box, mtot is the total mass, and g is an exponent that is intended to give more weight to 
either high density (when g > 0) or low density regions (g < 0). To improve the statistics, 
the resnlt is averaged over 1024 realizations obtained by shifting the origin of the system 
by mnltiples of the grid spacing and taking into acconnt the periodicity of the bonndary 
conditions. 


A few examples of compntation of Dq are shown in Fig. In practice, Dq is given by 
the slope of the cnrves for intermediate length scales. Note that for large ln(£) the slope is 
always eqnal to nnity, signalling that the system becomes homogeneons at large scales. The 
nncertainty in the linear regression procednre yields the error bars that appear in Fig. 

It can be proven [33] that Dq shonld be a monotonically decreasing (or fiat) fnnction of 
the exponent g. However, N-body simulations showed that, while Dq displays the expected 
trend for g > 0, it is an increasing function for g < 0 (open circles in Fig. [^. Since negative 
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no cut q= -3.0 
no cut q= -1.0 
no cut q= 3.0 
cut 10“^ q= -3.0 
cut 10“^ q=-1.0 
cut 10“^ q= 3.0 
cut 10“^ q= -3.0 
cut 10“^ q=-1.0 
cut 10“^ q= 3.0 
cut 1 q= -3.0 
cut 1 q= -1.0 
cut 1 q= 3.0 


FIG. 6: (Color online). Computation of the fractal dimension Dq from Eq. (25), for three values 
of q = -3,-1 and 3, and several cut-offs; pth = 10“^, 10“^, and 1. The thick lines indicate the 
ranges over which the slope was computed (these ranges are the same for different values of g at a 
given cut-off). 


values of q overrepresent low-density regions, this behavior was attributed to poor sampling 
of these regions, where the number of particles is small and the statistics noisy. 

Vlasov codes, by sampling the entire phase space with the same accuracy irrespective of the 
matter content, should provide better results precisely in such low-density regions. This is 
indeed what we observe on Fig. For positive values of q, which are dominated by large- 
density regions, the Vlasov and N-body results are in agreement; in contrast, for negative 
q the Vlasov results (open squares) level off at Dq ^ 1 |3H]. At face value, these findings 
suggest that the matter distribution is fractal at high densities (because Dq < 1 for q > 0), 
whereas it is homogeneous at low densities {Dq ^ 1 for g < 0). If confirmed, this would 
be an important result for our understanding of the distribution of matter in the universe. 
Nevertheless, it cannot be ruled out that the leveling off of Dq is partly due to numerical 
diffusion, which washes out the small-scale structures that develop over time. 

To understand why N-body codes fail to reproduce correctly the g < 0 region, we introduced 
an artificial cut-off in the density p issued from the Vlasov simulations. Thus, values of p that 
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FIG. 7: (Color online). Fractal dimension Dq for various values of the exponent q and different cut¬ 
offs of the matter density. Open circles correspond to the N-body results obtained with N ~ 262 000 
particles. The black dashed curve corresponds to the full Vlasov results (no cut-off). Other lines 
correspond to Vlasov results with cut-off threshold at pmin = 10“'^, 10“^, 1, and 4. 

are below a certain threshold are set to zero. This is the same procedure that was applied 
earlier to the power spectrum. In Fig. we show the results for four values of the threshold, 
Pmin = 10“"^, 10“^, 1, and 4 (note that, although p is normalized to unity, the fluctuations 
can be much larger, as seen in Fig. [^. It is clear that, by increasing the cut-off level, the 
Vlasov results progressively move towards the N-body results. Interestingly, we observe that 
the Vlasov and N-body results start to coincide for a threshold value pth ~ 4. This is in 
agreement with the cut-off value of the density that is necessary to recover the slope of the 
power spectrum observed in the N-body code (see Fig. and related discussion). These 
hndings strongly suggest that the incorrect behavior of Dq observed in N-body simulations 
is indeed due to poor sampling of the low-density regions, and that this drawback can be 
overcome by using a numerical approach based on a uniform meshing of the phase space. 
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V. DISCUSSION AND CONCLUSIONS 


Most numerical simulations of self-gravitating systems are performed using N-body codes, 
which solve the equations of motion of a large number of interacting particles. This is an 
approximation, since the number of “particles” in a real system is virtually infinite, whereas 
simulations are limited to a few million particles. Ideally, one should instead solve the 
Vlasov-Poisson equations, but this is more costly in terms of memory storage and computing 
time. Except for some simplified cases [SU |35] , Vlasov simulations are out of reach of current 
computer capabilities for 3D problems, although they are now feasible for ID problems. Here, 
we have shown an application of Vlasov codes to cosmological simulations of an expanding 
universe. A key point was the choice of the most suitable scaling factors, which map the 
original phase space (a:, v) onto a scaled phase space (^, rj) that is bounded for all times, 
thus optimizing the number of mesh points. 

The results confirmed the appearance of self-similar clustering in the phase space and a 
power-law spectrum similar to that observed for N-body simulations. The box-counting 
fractal dimension Dg is flat and close to unity for g < 0 and decreasing for g > 0, suggesting 
that the matter distribution is not the same at low and high densities. This different behavior 
was also visible in the power spectra observed at low and high densities, and may offer an 
insight into our understanding of large cosmological structures such as voids. Nevertheless, 
these preliminary results, which are potentially sensitive to the numerical resolution of the 
code, would require more studies to be fully confirmed. Other approaches based on equal 
mass partitions [36] may provide further information on the low-density regions. 


Appendix A: New scaling 

It was observed in N-body numerical simulations that the variance of the scaled velocity 
(“thermal” velocity) rjth = grows exponentially in time. This is a problem for grid- 

based Vlasov simulations, since one would need to mesh a very large velocity space in order to 
keep the distribution function inside the computational box for all times. Therefore, we want 
to look for a new scaling for which the scaled velocity is bound in time, i.e., (? 7 ^) ~ const. 
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More precisely, N-body numerical simulations show that (see Fig. [^: 


r]th = Vo exp ( -ujj Ooid 


(Al) 


where 9oid is the scaled time obtained with the standard scaling (/3 = 1). Using the time t 
and remembering that aujto = 1 , we obtain 


t \ 

Vth = Vo \ - ] 


(A 2 ) 


The relationship between the thermal velocities Vth and Vth is deduced from Eq. (|^, where 
we neglect the last term because C decreases to zero with time: Vth = {C/A^)vth- Therefore, 


using Eq. (A2), we obtain: 




Uh ^2 Vth 


-1/3 


Vo 


^ \ l/3o f t\ 


= Vo 


to 


(A3) 


Ao J \to J 

where we have used the standard scaling exponents (3 = 1 and 7 = 2/3. 

Now we want to hnd a new scaling where Vth is bounded in time. In Eq. (|^, we still keep 
the exponent 7 = 2/3 (because it represents the physical expansion rate of an Einstein-de 
Sitter universe), and look for an exponent (3 that yields Vth ^ const. We have: 


A 2 


Vth = ^ vth = \ — 


C 


^ N /3-(2/3) / ^ N il-a)/3a 


to 




In order for Vth to be constant in time, one needs to satisfy: 


which yields 


/5 = 


3a 3 


3a - 1 9 -V 2 


3a 


9 


0.84. 


(A4) 


(AS) 


(A 6 ) 


This is the value of (3 used in our simulations. 
The new scaling derived above 

9-\/2 


, N 2/3 


(A7) 


yields a different equation of motion, where some of the coefficients are time-dependent. 
Substituting Eq. (A7) into the general form of the equation of motion ([^, one obtains 


d^^ 3fo Vto/ d6' 


9tl 


to) 


2(/3-l) / . \ 2(/3-l) 


(AS) 
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FIG. 8: Evolution of the scaled thermal velocity r]th as a function of the scaled time Ooid normalized 
to (standard scaling with j3 = 1) for an N-body simulation with initial condition similar to 
that used for the Vlasov case. The measured slope is 0.361, close to the theoretical value 1/3 
(straight line). 


where now jS = ^ ^ . The relationship between the times t and 6 is obtained by integrating 


Eq. (|^ with the condition that 9 = Q when t = This yields 


t 

to 


ti 


o \ l/(l-/3) 

-0 + 1 




(A9) 


Substituting the above expression into Eq. (A8), we obtain 
d^^ 4-3/3 1 


d02 


+ 


1 d^ 2 1 

(l-/3)^ + ld0“^[(l-/3)|- + l ]2 




[(l-/^)g + iP 


S 


(AlO) 


Using the relation (A6) and remembering that auji^ = 1, one gets: 
d^^ a + 1 


d02 


+ 




de 


2, ,2 


a ui 


e = 


£ 


3 u;j,9/3 + 1 d0 9 [a;j0/3 + 1]2 ^ [u;j9/3 + 1]2 ’ 


Dehning £ = £/oj‘] and 9 = ojj9, Eq. (All) becomes 


de 


a 


•e = 


£ 


d ^ a + 1 1 

d 02 3 9/3 + 1 S 9 [e/3 + 1]2 ^ ^ [0/3 + 1]2 ’ 


(All) 


(A12) 
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and finally, with a = 3/v^, 

3 + ^2 


de 


d02 3^2 0/3 + 1 d0 [0/3 + 1]2 

which is identical to Eq. (8) in the main text. 




E 


[0/3 + l]2’ 


(A13) 


Appendix B: Relation between the old and new scaled times 


Let us call Odd the scaled time for which A^{t) = t/to (standard scaling, (3 = 1), and 6new 
the scaled time for which A{t)‘^ = with (3 = (new scaling). For the standard 

and new scalings, the relationships between the scaled times Ooid and 9new and the real time 
t read as follows: 

i (9oid 

- = exp — 

^0 \ ^0 


and 


t 

to 


l + (l-/9) 


6r 


to 


l/(l-/3) 


We take the same to in both cases, since it is the instant at which the real and scaled times 
coincide. Equating the two expressions for t/to, we hnd: 




(Bl) 


or, normalizing the scaled times to uj: 

^j9oid = 


In [1 + q;(1 — (3)uj9 

new\ 1 


(B2) 


a{l — j3) 

with a{l — f3) = 1/3. 

For instance, when u)j9new = 300, we obtain uj9oid ~ 13.85. This value is in accordance 
with the simulation times used for the N-body simulations using the “old” scaling variables. 
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